rm(list = ls())
library(pacman)
p_load(lme4,lmerTest ,performance, tidyverse,lubridate, lattice, devtools, VGAM, ggpubr, plotly, knitr)
rm(list = ls())
 setwd("F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/dataset/data_sep")

load(file="processed_data.RData")

# check the Missing data on test.1000 
miss.test <- data0 %>% filter(is.na(test.thousand) ==T) %>% 
  select(iso_code, country, date.date, test.thousand ) 
miss.test.country <- data0 %>% filter(is.na(test.thousand) ==T) %>% 
  distinct(iso_code)
## Countries dont have the missing data
no.mist.test <- data0 %>% filter(is.na(test.thousand) ==F) %>% 
  distinct(iso_code, country)  # 31 countries
save(no.mist.test, file = "no.miss.test.RData")

# Check missing of week.case.eff 
miss.case <- data0 %>% filter(is.na(total.cases) == T)%>% 
  select(iso_code, country, date.date, week.case.eff )
miss.case.cnt <- data0 %>% filter(is.na(total.cases) == T)%>% 
  distinct(iso_code)
setwd("F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/dataset/data_sep")
vac.week20 <- data0.anal %>% filter(time == 19) %>% select(iso_code, country, total.cases,tt.vac.100, pp.vac.100, ful.vac.100)
write.table(vac.week20, file="vac_week20_plotdata.csv", sep = ",")

save(vac.week20, file = "vac_week20.RData")
# check.data <- data0.anal %>% select(id, iso_code, time, date.date, tt.vac.100, pp.vac.100, vac.100.lag, ful.vac.100,  test.thousand, total.cases,  week.case.eff) 
# 
# 
# check.data <- data0.ful %>% select(id, iso_code, time, date.date, tt.vac.100, pp.vac.100, ful.vac.100,  test.thousand, avg.string, week.case.eff) 
# 
# check.data <- vac.period %>% select(id, iso_code, date.date, population, tt.vac.100, pp.vac.100, ful.vac.100, total.cases, test.thousand) 
# 
# check.vac <- data0.anal %>% filter(is.na(vac.100.cate) == T) %>% select(iso_code, id, time, pp.vac.100, vac.100.cate)
# hist(data0.anal$week.case.eff
#      )
# check <- data0.anal %>% filter(week.case.eff <0) %>% select(iso_code, id, week.case.eff)
# check <- data0.anal %>% filter(iso_code == "KAZ") %>% 
#   select(iso_code, date.date, time, total.cases, week.case.eff
#          )


# KAZ data wrong
check.time <- data0.anal %>% select(id, week.case.eff, growth.rate.manual, timeStart, time)
## Adding missing grouping variables: `iso_code`
 # Group plot
plot.group.case <- ggplot(data0.anal, aes(time+2, week.case.eff)) +
  stat_smooth(method = "loess", se = F, size = 2) +
  stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
  xlab("Week") +
  ylab(" Average weekly growth rate")+ 
  theme_bw() + 
  theme(strip.background = element_blank()) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(legend.position = "bottom") + 
  scale_x_continuous(breaks = seq(2, 25, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.case
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing non-finite values (stat_summary).

plot.group.string <- ggplot(data0.anal, aes(time, string.idx)) +
  stat_smooth(method = "loess", se = F, size = 2) +
  stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
  xlab("Week") +
  ylab(" Stringency index")+ 
  theme_bw() + 
  theme(strip.background = element_blank()) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 23, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.string 
## `geom_smooth()` using formula 'y ~ x'

cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig1_NPI_vac.period_vac.period_Sep27.pdf",
          width = 18, height = 8, #inch
          onefile = TRUE, family = "Segoe UI")

merge1 <- ggarrange(plot.group.string, plot.group.case, hjust = -1,
          ncol = 2, nrow = 1)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).

## Warning: Removed 3 rows containing non-finite values (stat_summary).
merge1

dev.off()
## png 
##   2
ggplotly(plot.group.case)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).

## Warning: Removed 3 rows containing non-finite values (stat_summary).
ggplotly(plot.group.string)
## `geom_smooth()` using formula 'y ~ x'
#https://www.guru99.com/r-scatter-plot-ggplot2.html#2
# clean daily grow rate = 0, and NAs
## NAs = 36 obs 
## zero growth rate = 58

## Check zero growth rate
data0.anal1 <- data0.anal
check <- data0.anal1 %>% filter(week.case.eff == 0 )  
check1 <- data0.anal1 %>% filter(week.case.eff != 0 ) %>% select(id, week.case.eff) 
## Adding missing grouping variables: `iso_code`
## Check NAs
summary(data0.anal$week.case.eff)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## 0.000000 0.002658 0.004951 0.010112 0.009787 0.357987        3
#probit transformation
data0.anal <- data0.anal %>% filter(week.case.eff != 0) %>% mutate(probit.case = probitlink( week.case.eff))

qqnorm(data0.anal$probit.case)

hist(data0.anal$probit.case)

boxplot(data0.anal$probit.case)

summary(data0.anal$test.policy.red)
## Warning: Unknown or uninitialised column: `test.policy.red`.
## Length  Class   Mode 
##      0   NULL   NULL
data0.anal <- data0.anal %>% mutate(test.policy.red = case_when(
  test.policy==0~0,
  test.policy==1~0,
  test.policy==2~0, 
  test.policy==3~1
))
data0.anal$test.policy.red <- factor(data0.anal$test.policy.red)



data0.anal$pop_log <- log(data0.anal$population)
hist(data0.anal$pop_log)

data0.anal$pop_denlog <- log(data0.anal$pop_density)
hist(data0.anal$pop_denlog)

data0.anal$median_age.log <- log(data0.anal$median_age)
hist(data0.anal$median_age.log)

data0.anal$gdp.log <- log(data0.anal$gdp)
hist(data0.anal$gdp.log)

data0.anal$incomegroup1 <- factor(data0.anal$incomegroup, levels = c("Low income", "Lower middle income", "Upper middle income" , "High income"))
data0.anal <- data0.anal %>% mutate(incomegroup2 = case_when(
  incomegroup1 == "Low income" ~ "LMIC",
  incomegroup1 == "Lower middle income" ~ "LMIC",
  incomegroup1 == "Upper middle income" ~ "upper_income",
  incomegroup1 == "High income" ~ "high_income"
)) 
data0.anal$incomegroup2 <- factor(data0.anal$incomegroup2, levels = c( "LMIC", "upper_income","high_income"))

# 
# Sig at >30
# data0.anal <- data0.anal %>% mutate(vac.bin = case_when(
#   pp.vac.100 <=5 ~ 0, 
#   pp.vac.100 <=10 & pp.vac.100 >5 ~ 1,
#   pp.vac.100 <=30 & pp.vac.100 >10 ~ 2,
#   pp.vac.100 > 30 ~ 3,
# ))
# data0.anal <- data0.anal %>% mutate(vac.bin = case_when(
#   pp.vac.100 <=1 ~ 0,
#   pp.vac.100 <=8 & pp.vac.100 >1 ~ 1,
#   pp.vac.100 <=16 & pp.vac.100 >8 ~ 2,
#   pp.vac.100 > 16 ~ 3,
# ))

data0.anal$pp.vac.1 <- data0.anal$pp.vac.100

data0.anal <- data0.anal %>% mutate(pp.vac.100 = case_when(
  pp.vac.1 <1 ~ 0,
  pp.vac.1 <5 & pp.vac.1 >=1 ~ 1,
  pp.vac.1 <10 & pp.vac.1 >=5 ~ 2,
  pp.vac.1 <30 & pp.vac.1 >=10 ~ 3,
  pp.vac.1 > 30 ~ 4,
))



data0.anal$pp.vac.100 <- factor(data0.anal$pp.vac.100)

summary(data0.anal$ful.vac.100)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.670   4.685   3.185  60.780
data0.anal <- data0.anal %>% mutate(ful.100.cate = case_when(
  ful.vac.100 <=10 & ful.vac.100 >=0 ~ 0, 
   ful.vac.100 <=30 & ful.vac.100 >10~ 1,
  ful.vac.100 <100 & ful.vac.100 >30 ~ 2
))
data0.anal$ful.100.cate <- factor(data0.anal$ful.100.cate)
library(performance)
aic <- function(model){
  fix2 <- performance(model)
  x<- as.character(substitute(model))
  fix2 <- fix2 %>% 
    mutate( npi = x)
  }
model.school.cls <- lmer(probit.case ~ school.close.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.school.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ school.close.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     87.5    122.1    -35.8     71.5      547 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7392 -0.5252 -0.0208  0.4549  7.4866 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2254225 0.47479       
##           time        0.0008684 0.02947  -0.63
##  Residual             0.0477187 0.21845       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       -2.613e+00  9.421e-02  3.133e+01 -27.734  < 2e-16 ***
## school.close.red2 -2.283e-02  4.198e-02  5.399e+02  -0.544  0.58687    
## school.close.red3  1.118e-01  3.445e-02  5.414e+02   3.246  0.00124 ** 
## time              -7.892e-04  5.806e-03  2.765e+01  -0.136  0.89286    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sch..2 sch..3
## schl.cls.r2 -0.169              
## schl.cls.r3 -0.233  0.549       
## time        -0.615 -0.035 -0.002
x1 <- aic(model.school.cls)


model.work.cls <- lmer(probit.case ~ work.close.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.work.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ work.close.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    104.8    139.3    -44.4     88.8      547 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3332 -0.5075 -0.0527  0.5271  7.2350 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2281360 0.47764       
##           time        0.0009398 0.03066  -0.63
##  Residual             0.0491885 0.22178       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      -2.575427   0.096201  33.054626 -26.771   <2e-16 ***
## work.close.red2   0.034585   0.043643 549.693377   0.792    0.428    
## work.close.red3   0.025341   0.047098 547.950575   0.538    0.591    
## time             -0.001756   0.006040  28.153571  -0.291    0.773    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) wrk..2 wrk..3
## wrk.cls.rd2 -0.271              
## wrk.cls.rd3 -0.257  0.670       
## time        -0.594 -0.072 -0.040
x2<- aic(model.work.cls)



model.pubevent <- lmer(probit.case ~ cancel.pubevent.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pubevent)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ cancel.pubevent.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    103.3    133.6    -44.7     89.3      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3111 -0.5075 -0.0424  0.5259  7.1903 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2274454 0.47691       
##           time        0.0009387 0.03064  -0.64
##  Residual             0.0492832 0.22200       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -2.565247   0.098900  36.761623 -25.938   <2e-16 ***
## cancel.pubevent.red2   0.012776   0.043990 546.532929   0.290    0.772    
## time                  -0.001423   0.006021  27.525849  -0.236    0.815    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnc..2
## cncl.pbvn.2 -0.368       
## time        -0.599 -0.006
x3<- aic(model.pubevent)


model.restr.gather <- lmer(probit.case ~ restric.gather.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.gather)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ restric.gather.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     94.6    133.5    -38.3     76.6      546 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1704 -0.4927 -0.0354  0.5313  6.7138 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  iso_code (Intercept) 0.235730 0.48552       
##           time        0.001022 0.03197  -0.65
##  Residual             0.047838 0.21872       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          -2.302641   0.119205  70.361538 -19.317  < 2e-16 ***
## restric.gather.red1  -0.296975   0.110053 542.504290  -2.698 0.007183 ** 
## restric.gather.red3  -0.218769   0.080203 523.422597  -2.728 0.006592 ** 
## restric.gather.red4  -0.276462   0.079468 521.500973  -3.479 0.000546 ***
## time                 -0.000747   0.006262  27.758479  -0.119 0.905896    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rst..1 rst..3 rst..4
## rstrc.gth.1 -0.464                     
## rstrc.gth.3 -0.576  0.686              
## rstrc.gth.4 -0.613  0.674  0.878       
## time        -0.497 -0.036 -0.035 -0.031
x4 <- aic(model.restr.gather)



model.cls.trans <- lmer(probit.case ~ close.transport.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.cls.trans)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ close.transport.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     86.7    116.9    -36.3     72.7      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1660 -0.5231 -0.0299  0.5296  6.9448 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  iso_code (Intercept) 0.247529 0.4975        
##           time        0.001024 0.0320   -0.65
##  Residual             0.047302 0.2175        
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -2.450992   0.098912  31.484667 -24.780  < 2e-16 ***
## close.transport.red1  -0.174516   0.042036 549.166296  -4.152 3.83e-05 ***
## time                  -0.000759   0.006262  27.902486  -0.121    0.904    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cls..1
## cls.trnsp.1 -0.252       
## time        -0.628 -0.025
x5 <- aic(model.cls.trans)


model.stay.home <- lmer(probit.case ~ stay.home.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.stay.home)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ stay.home.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    103.3    133.5    -44.6     89.3      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2604 -0.5068 -0.0481  0.5328  7.1337 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2295284 0.47909       
##           time        0.0009564 0.03093  -0.64
##  Residual             0.0492192 0.22185       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     -2.547814   0.094101  30.050805 -27.075   <2e-16 ***
## stay.home.red1  -0.013358   0.035099 546.888338  -0.381    0.704    
## time            -0.001325   0.006077  27.999834  -0.218    0.829    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sty..1
## stay.hm.rd1 -0.192       
## time        -0.629 -0.038
x6 <- aic(model.stay.home)


model.restr.move <- lmer(probit.case ~ restric.move.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.move)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ restric.move.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    101.8    132.0    -43.9     87.8      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3178 -0.5154 -0.0484  0.5343  7.1029 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2400995 0.49000       
##           time        0.0009943 0.03153  -0.65
##  Residual             0.0489042 0.22114       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -2.534366   0.095680  28.964016 -26.488   <2e-16 ***
## restric.move.red2  -0.046292   0.036008 547.304999  -1.286    0.199    
## time               -0.001214   0.006183  27.635482  -0.196    0.846    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rst..2
## rstrc.mv.r2 -0.165       
## time        -0.643 -0.025
x7 <- aic(model.restr.move)


model.inter.trav <- lmer(probit.case ~ inter.travel +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.inter.trav)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ inter.travel + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    104.3    143.2    -43.2     86.3      546 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3336 -0.4979 -0.0561  0.5121  7.1584 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2363377 0.48615       
##           time        0.0009839 0.03137  -0.65
##  Residual             0.0488493 0.22102       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    -2.474062   0.109431  49.579318 -22.608   <2e-16 ***
## inter.travel2  -0.106001   0.061352 547.033631  -1.728   0.0846 .  
## inter.travel3  -0.072392   0.068242 553.914190  -1.061   0.2892    
## inter.travel4  -0.088579   0.077671 550.851424  -1.140   0.2546    
## time           -0.001490   0.006165  27.593026  -0.242   0.8109    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) intr.2 intr.3 intr.4
## inter.trvl2 -0.465                     
## inter.trvl3 -0.475  0.742              
## inter.trvl4 -0.457  0.717  0.696       
## time        -0.571  0.006 -0.003  0.044
x8 <- aic(model.inter.trav)

model.test <- lmer(probit.case ~ test.policy.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ test.policy.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.2    132.4    -44.1     88.2      546 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2927 -0.5068 -0.0320  0.5216  7.1756 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2306891 0.48030       
##           time        0.0009633 0.03104  -0.64
##  Residual             0.0490835 0.22155       
## Number of obs: 553, groups:  iso_code, 28
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       -2.562544   0.096770  32.166052 -26.481   <2e-16 ***
## test.policy.red1   0.014886   0.050557 525.740203   0.294    0.769    
## time              -0.001444   0.006101  27.889687  -0.237    0.815    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tst..1
## tst.plcy.r1 -0.291       
## time        -0.608 -0.055
x9 <- aic(model.test)

model.tracing <- lmer(probit.case ~ contact.tracing +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.tracing)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ contact.tracing + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    105.0    139.4    -44.5     89.0      542 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3282 -0.5107 -0.0292  0.5465  7.2002 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2285061 0.47802       
##           time        0.0009832 0.03136  -0.64
##  Residual             0.0490618 0.22150       
## Number of obs: 550, groups:  iso_code, 28
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      -2.647e+00  1.439e-01  1.302e+02 -18.388   <2e-16 ***
## contact.tracing1  7.209e-02  1.133e-01  5.272e+02   0.636    0.525    
## contact.tracing2  1.079e-01  1.159e-01  5.186e+02   0.931    0.352    
## time             -8.887e-04  6.165e-03  2.810e+01  -0.144    0.886    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cntc.1 cntc.2
## cntct.trcn1 -0.744              
## cntct.trcn2 -0.762  0.930       
## time        -0.458  0.047  0.062
x10 <- aic(model.tracing)



model.mask <- lmer(probit.case ~ wear.mask +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mask)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ wear.mask + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    109.1    148.0    -45.6     91.1      544 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3215 -0.5246 -0.0494  0.5305  7.1577 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2281657 0.47767       
##           time        0.0009572 0.03094  -0.64
##  Residual             0.0494212 0.22231       
## Number of obs: 553, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -2.621808   0.374853 29.197539  -6.994 1.05e-07 ***
## wear.mask2   0.073485   0.394960 33.419385   0.186    0.854    
## wear.mask3   0.081558   0.377430 27.963831   0.216    0.830    
## wear.mask4   0.063663   0.377380 27.949729   0.169    0.867    
## time        -0.001686   0.006118 28.510302  -0.276    0.785    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) wr.ms2 wr.ms3 wr.ms4
## wear.mask2 -0.927                     
## wear.mask3 -0.966  0.949              
## wear.mask4 -0.968  0.952  0.993       
## time       -0.161  0.017 -0.007  0.006
x11 <- aic(model.mask)

model.info <- lmer(probit.case ~ info.campaign +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.info)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ info.campaign + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     96.6    131.2    -40.3     80.6      547 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0543 -0.5110 -0.0343  0.5388  7.0399 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  iso_code (Intercept) 0.234285 0.48403       
##           time        0.001039 0.03224  -0.69
##  Residual             0.048432 0.22007       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     -2.665876   0.251606 465.932475 -10.595   <2e-16 ***
## info.campaign1   0.387860   0.251629 525.319857   1.541    0.124    
## info.campaign2   0.105837   0.233361 506.472699   0.454    0.650    
## time            -0.002833   0.006328  27.825543  -0.448    0.658    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) inf.c1 inf.c2
## info.cmpgn1 -0.863              
## info.cmpgn2 -0.929  0.927       
## time        -0.270 -0.015  0.015
x12 <- aic(model.info)
#population: no
model.pop <- lmer(probit.case ~ pop_log +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pop)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ pop_log + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    103.3    133.5    -44.6     89.3      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2884 -0.5095 -0.0432  0.5308  7.1668 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2202337 0.46929       
##           time        0.0009495 0.03081  -0.62
##  Residual             0.0492504 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept) -2.264569   0.679272 28.403577  -3.334  0.00239 **
## pop_log     -0.017000   0.039448 28.004982  -0.431  0.66980   
## time        -0.001411   0.006052 28.044731  -0.233  0.81741   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) pop_lg
## pop_log -0.991       
## time    -0.084  0.000
y1 <- aic(model.pop)

#population density

model.popden <- lmer(probit.case ~ pop_denlog +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.popden)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ pop_denlog + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     97.6    127.8    -41.8     83.6      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2930 -0.5117 -0.0337  0.5312  7.1691 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.1776080 0.42144       
##           time        0.0009495 0.03081  -0.62
##  Residual             0.0492504 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -2.008411   0.223225 30.842406  -8.997 3.93e-10 ***
## pop_denlog  -0.107027   0.040693 27.997629  -2.630   0.0137 *  
## time        -0.001406   0.006053 28.044514  -0.232   0.8180    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) pp_dnl
## pop_denlog -0.931       
## time       -0.231  0.000
y2 <- aic(model.popden)

# median age
model.age <- lmer(probit.case ~ median_age +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ median_age + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.2    132.5    -44.1     88.2      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2915 -0.5070 -0.0349  0.5270  7.1630 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2162986 0.46508       
##           time        0.0009495 0.03081  -0.63
##  Residual             0.0492502 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -2.201318   0.333172 29.660521  -6.607 2.74e-07 ***
## median_age  -0.011469   0.010413 28.012157  -1.101    0.280    
## time        -0.001411   0.006052 28.044495  -0.233    0.817    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) medn_g
## median_age -0.963       
## time       -0.173  0.000
y3 <- aic(model.age)


# gdp
model.gdp <- lmer(probit.case ~ gdp.log +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.gdp)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ gdp.log + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    100.9    131.2    -43.5     86.9      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2935 -0.5043 -0.0357  0.5266  7.1522 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2319973 0.48166       
##           time        0.0009494 0.03081  -0.68
##  Residual             0.0492505 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept) -1.526549   0.631479 28.576159  -2.417   0.0222 *
## gdp.log     -0.106908   0.064952 27.999482  -1.646   0.1110  
## time        -0.001418   0.006052 28.045002  -0.234   0.8164  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) gdp.lg
## gdp.log -0.989       
## time    -0.102  0.000
y4 <- aic(model.gdp)

# test per 1000 => too many missing because the early stage of pandemic
model.test1000 <- lmer(probit.case ~ test.thousand +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test1000)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ test.thousand + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.3    132.6    -44.2     88.3      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2922 -0.5113 -0.0432  0.5331  7.1611 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  iso_code (Intercept) 0.233387 0.48310       
##           time        0.000932 0.03053  -0.66
##  Residual             0.049246 0.22191       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   -2.525e+00  9.730e-02  3.182e+01 -25.949   <2e-16 ***
## test.thousand -9.077e-05  8.569e-05  2.966e+01  -1.059    0.298    
## time          -6.059e-05  6.135e-03  3.011e+01  -0.010    0.992    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tst.th
## test.thosnd -0.290       
## time        -0.561 -0.208
y5 <- aic(model.test1000)


# percentage of urban
model.per_urban <- lmer(probit.case ~ per_urban + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.per_urban)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ per_urban + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    100.0    130.2    -43.0     86.0      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2963 -0.5018 -0.0356  0.5263  7.1450 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2407254 0.49064       
##           time        0.0009493 0.03081  -0.71
##  Residual             0.0492508 0.22193       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -2.254156   0.177031 35.178733 -12.733 9.86e-15 ***
## per_urban   -0.004980   0.002481 28.003151  -2.007   0.0545 .  
## time        -0.001422   0.006052 28.044844  -0.235   0.8159    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) pr_rbn
## per_urban -0.846       
## time      -0.382  0.000
y6 <- aic(model.per_urban)

# uhc_index
model.uhc <- lmer(probit.case ~ uhc_index + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.uhc)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ uhc_index + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    103.0    133.2    -44.5     89.0      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2921 -0.5068 -0.0358  0.5279  7.1601 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2253242 0.47468       
##           time        0.0009495 0.03081  -0.64
##  Residual             0.0492503 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -2.297240   0.408576 29.177211  -5.623 4.41e-06 ***
## uhc_index   -0.003871   0.005987 28.004973  -0.647    0.523    
## time        -0.001413   0.006052 28.044511  -0.233    0.817    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) uhc_nd
## uhc_index -0.975       
## time      -0.145  0.000
y7 <- aic(model.uhc)

# health worker density
model.hworker <- lmer(probit.case ~ healthworker_den + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.hworker)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ healthworker_den + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.0    132.2    -44.0     88.0      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2972 -0.5022 -0.0356  0.5291  7.1524 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2270257 0.47647       
##           time        0.0009495 0.03081  -0.66
##  Residual             0.0492503 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      -2.404549   0.153802 35.176866 -15.634   <2e-16 ***
## healthworker_den -0.002500   0.002054 28.006164  -1.217    0.234    
## time             -0.001415   0.006052 28.044592  -0.234    0.817    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hlthw_
## hlthwrkr_dn -0.802       
## time        -0.398  0.000
y8 <- aic(model.hworker)

#Income group
model.income <- lmer(probit.case ~ incomegroup2 + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.income)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ incomegroup2 + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.0    136.6    -43.0     86.0      547 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2934 -0.5038 -0.0287  0.5261  7.1545 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2221627 0.47134       
##           time        0.0009495 0.03081  -0.68
##  Residual             0.0492501 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)              -2.443163   0.115400 35.758920 -21.171   <2e-16 ***
## incomegroup2upper_income -0.102084   0.164078 28.022348  -0.622    0.539    
## incomegroup2high_income  -0.300860   0.157239 27.999869  -1.913    0.066 .  
## time                     -0.001415   0.006052 28.044568  -0.234    0.817    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) incmgrp2p_ incmgrp2h_
## incmgrp2pp_ -0.498                      
## incmgrp2hg_ -0.519  0.365               
## time        -0.541  0.001      0.000
y9 <- aic(model.income)

# SDI 
model.sdi <- lmer(probit.case ~ sdi_index + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.sdi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ sdi_index + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    101.7    131.9    -43.8     87.7      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2950 -0.5065 -0.0335  0.5297  7.1547 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2254544 0.47482       
##           time        0.0009495 0.03081  -0.66
##  Residual             0.0492504 0.22192       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -2.140177   0.324088 29.971522  -6.604 2.63e-07 ***
## sdi_index   -0.626891   0.470193 28.002484  -1.333    0.193    
## time        -0.001416   0.006052 28.044655  -0.234    0.817    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) sd_ndx
## sdi_index -0.959       
## time      -0.189  0.000
y10 <- aic(model.sdi)

# mobi.composite 
model.mobi <- lmer(probit.case ~ mobi.composite + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.mobi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ mobi.composite + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    -39.8    -10.9     26.9    -53.8      451 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2721 -0.5685 -0.0579  0.4959  3.5716 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  iso_code (Intercept) 0.140904 0.37537       
##           time        0.001008 0.03175  -0.57
##  Residual             0.036810 0.19186       
## Number of obs: 458, groups:  iso_code, 23
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    -2.697e+00  8.170e-02  2.471e+01 -33.011   <2e-16 ***
## mobi.composite  7.427e-04  1.003e-03  4.491e+02   0.741    0.459    
## time            5.051e-03  6.813e-03  2.179e+01   0.741    0.466    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mb.cmp
## mobi.compst  0.193       
## time        -0.562  0.056
y11 <- aic(model.mobi)

# pp.vac.100


vac.100 <- lmer(probit.case ~ pp.vac.100 + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(vac.100)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ pp.vac.100 + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.4    145.6    -41.2     82.4      545 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2500 -0.4979 -0.0390  0.5044  7.0300 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2418162 0.49175       
##           time        0.0009176 0.03029  -0.65
##  Residual             0.0485558 0.22035       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  -2.531025   0.096004  29.241183 -26.364  < 2e-16 ***
## pp.vac.1001  -0.060768   0.040118 529.889865  -1.515  0.13043    
## pp.vac.1002  -0.117374   0.058890 549.454377  -1.993  0.04675 *  
## pp.vac.1003  -0.163785   0.074183 553.142389  -2.208  0.02766 *  
## pp.vac.1004  -0.256442   0.097664 506.277510  -2.626  0.00891 ** 
## time          0.006671   0.006734  44.127606   0.991  0.32725    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p..1001 p..1002 p..1003 p..1004
## pp.vac.1001 -0.161                                
## pp.vac.1002 -0.135  0.724                         
## pp.vac.1003 -0.108  0.684   0.806                 
## pp.vac.1004 -0.087  0.556   0.656   0.745         
## time        -0.520 -0.313  -0.383  -0.442  -0.429
y11 <- aic(vac.100)

# ful.vac.100


ful.vac.100 <- lmer(probit.case ~ ful.vac.100 + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(ful.vac.100)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ ful.vac.100 + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     91.4    121.7    -38.7     77.4      548 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3804 -0.4942 -0.0434  0.5060  7.2444 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2203538 0.46942       
##           time        0.0007844 0.02801  -0.60
##  Residual             0.0484613 0.22014       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  -2.576378   0.090747  28.201276 -28.391  < 2e-16 ***
## ful.vac.100  -0.007512   0.002124 398.790347  -3.537 0.000452 ***
## time          0.004567   0.005792  32.010772   0.788 0.436236    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) f..100
## ful.vac.100  0.068       
## time        -0.604 -0.292
y12 <- aic(ful.vac.100)

aic_table <- rbind(x1, x4, x5,x9 , y2, y6, y11, y12)
model.mul <- lmer(probit.case ~   close.transport.red + pp.vac.100 + school.close.red + pop_denlog + restric.gather.red + 
                    time  +  (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
                                            summary(model.mul)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ close.transport.red + pp.vac.100 + school.close.red +  
##     pop_denlog + restric.gather.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     49.6    123.0     -7.8     15.6      538 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7098 -0.5495  0.0116  0.4771  6.9168 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2146046 0.46325       
##           time        0.0009648 0.03106  -0.64
##  Residual             0.0426146 0.20643       
## Number of obs: 555, groups:  iso_code, 28
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -1.801946   0.252927  36.922632  -7.124 1.96e-08 ***
## close.transport.red1  -0.225375   0.049499 528.896720  -4.553 6.57e-06 ***
## pp.vac.1001           -0.090380   0.037939 527.396458  -2.382  0.01756 *  
## pp.vac.1002           -0.154076   0.055698 547.133298  -2.766  0.00586 ** 
## pp.vac.1003           -0.169693   0.070545 552.678279  -2.405  0.01648 *  
## pp.vac.1004           -0.291346   0.093020 525.947776  -3.132  0.00183 ** 
## school.close.red2      0.013337   0.041507 538.239817   0.321  0.74810    
## school.close.red3      0.184678   0.035316 535.362646   5.229 2.44e-07 ***
## pop_denlog            -0.102345   0.044287  28.754591  -2.311  0.02821 *  
## restric.gather.red1   -0.168341   0.118995 543.482660  -1.415  0.15773    
## restric.gather.red3   -0.089218   0.091998 548.902903  -0.970  0.33258    
## restric.gather.red4   -0.195709   0.089499 540.881089  -2.187  0.02919 *  
## time                   0.009285   0.006778  41.125757   1.370  0.17813    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
performance(model.mul)
## # Indices of model performance
## 
## AIC    |     BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
## ------------------------------------------------------------------
## 49.611 | 123.033 |      0.828 |      0.185 | 0.789 | 0.196 | 0.206
# #====================
 
# model.mul <- lmer(probit.case ~   close.transport.red + pp.vac.1 + school.close.red + pop_denlog + restric.gather.red + test.policy +
#                     time  +  (time|iso_code),
#                      data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
#                                             summary(model.mul)
# performance(model.mul1)
# 
# 
# 
# anova(model.mul, model.mul1) 
# # => No different select the simplifier model
#devtools::install_github("goodekat/redres")
library(redres)

data0.anal1 <- data0.anal %>%  filter( is.na(pp.vac.100) ==F) %>% filter(is.na(close.transport.red) == F) %>% filter(is.na(school.close.red) == F) %>%    filter(is.na(restric.gather.red) == F)  %>% 
  mutate(pop_denlog = log(pop_density))
 
rc_resids <- compute_redres(model.mul)
pm_resids <- compute_redres(model.mul, type = "pearson_mar")
## Warning in (subset(re_df, re_df$grp != "Residual")$sdcor)^2 * diag(q): longer
## object length is not a multiple of shorter object length
sc_resids <- compute_redres(model.mul, type = "std_cond")
## Warning in (subset(re_df, re_df$grp != "Residual")$sdcor)^2 * diag(q): longer
## object length is not a multiple of shorter object length
resids <- data.frame(data0.anal1$iso_code, rc_resids, pm_resids, sc_resids)
plot_redres(model.mul, type = "std_cond")
## Warning in (subset(re_df, re_df$grp != "Residual")$sdcor)^2 * diag(q): longer
## object length is not a multiple of shorter object length

plot_resqq(model.mul)

plot_ranef(model.mul)

plot(data0.anal1$time, rc_resids,ylim=c(-3,3))
abline(h=0, col="blue")

plot(data0.anal1$iso_code, sc_resids)
abline(h=c(0,2.96,-2.96), col="blue")

# # run other diagnosis
# p_load(patchwork, randomForest)
# 
 #check_model(model.mul) 
# 
 plot(check_distribution(model.mul))
## Warning: Removed 6 rows containing missing values (geom_segment).
## Warning: Removed 6 rows containing missing values (geom_point).

# 
 check_collinearity(model.mul)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                 Term  VIF Increased SE Tolerance
##  close.transport.red 1.54         1.24      0.65
##           pp.vac.100 1.34         1.16      0.74
##     school.close.red 1.22         1.10      0.82
##           pop_denlog 1.02         1.01      0.98
##   restric.gather.red 1.56         1.25      0.64
##                 time 1.25         1.12      0.80
# 
# pp_check(model.mul, 250)
library(glmmTMB)
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.3.4
## Current Matrix version is 1.3.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## 
## Attaching package: 'glmmTMB'
## The following object is masked from 'package:VGAM':
## 
##     betabinomial
# model.glm <- glmmTMB(week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red + pop_denlog + restric.gather.red + test.policy +
#                        time + (time|iso_code), data = data0.anal, family = beta_family(link="probit"), REML = TRUE)
# 
# summary(model.glm)

model.glm <- glmmTMB(week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red + pop_denlog + restric.gather.red +
                       time + (time|iso_code), data = data0.anal, family = beta_family(link="probit"), REML = TRUE)

summary(model.glm)
##  Family: beta  ( probit )
## Formula:          
## week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red +  
##     pop_denlog + restric.gather.red + time + (time | iso_code)
## Data: data0.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##  -4328.8  -4255.4   2181.4  -4362.8      551 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. Corr  
##  iso_code (Intercept) 0.1080878 0.32877        
##           time        0.0005248 0.02291  -0.60 
## Number of obs: 555, groups:  iso_code, 28
## 
## Dispersion parameter for beta family ():  234 
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.718147   0.200878  -8.553  < 2e-16 ***
## close.transport.red1 -0.150334   0.053208  -2.825 0.004722 ** 
## pp.vac.1001          -0.040151   0.039477  -1.017 0.309119    
## pp.vac.1002          -0.116616   0.054956  -2.122 0.033837 *  
## pp.vac.1003          -0.111118   0.068877  -1.613 0.106682    
## pp.vac.1004          -0.242739   0.095917  -2.531 0.011383 *  
## school.close.red2     0.001774   0.043782   0.041 0.967688    
## school.close.red3     0.128579   0.036852   3.489 0.000485 ***
## pop_denlog           -0.084909   0.035125  -2.417 0.015635 *  
## restric.gather.red1  -0.223756   0.125291  -1.786 0.074117 .  
## restric.gather.red3  -0.189786   0.098053  -1.936 0.052924 .  
## restric.gather.red4  -0.235337   0.092128  -2.554 0.010636 *  
## time                  0.006673   0.005550   1.202 0.229177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_load(plotly, modelbased)

data0.anal1$Predicted <- estimate_response(model.glm)$Predicted
## Warning in get_predicted.glmmTMB(model, data = data, predict = predict, : `predict = 'prediction'` is currently not available for glmmTMB models.
##   Changing to `predict = 'expectation'`.
plot2 <- data0.anal1 %>%
ggplot() +
geom_line(aes(x = log(week.case.eff), y = log(week.case.eff)), linetype = "dashed") +
geom_point(aes(x = log(week.case.eff) , y = log(Predicted), key=iso_code), color = "red") +
ylab("wADGR (predicted)") + xlab("wADGR (observed)")
## Warning: Ignoring unknown aesthetics: key
plot2

ggplotly(plot2, source = "select", tooltip = c("key") )
# Develop function marginal 
amex<- function(ame){
  x1 <- data.frame(list("AME" = ame*100, "factors" = as.character(substitute(ame))))
}
summary(model.glm)
##  Family: beta  ( probit )
## Formula:          
## week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red +  
##     pop_denlog + restric.gather.red + time + (time | iso_code)
## Data: data0.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##  -4328.8  -4255.4   2181.4  -4362.8      551 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. Corr  
##  iso_code (Intercept) 0.1080878 0.32877        
##           time        0.0005248 0.02291  -0.60 
## Number of obs: 555, groups:  iso_code, 28
## 
## Dispersion parameter for beta family ():  234 
## 
## Conditional model:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.718147   0.200878  -8.553  < 2e-16 ***
## close.transport.red1 -0.150334   0.053208  -2.825 0.004722 ** 
## pp.vac.1001          -0.040151   0.039477  -1.017 0.309119    
## pp.vac.1002          -0.116616   0.054956  -2.122 0.033837 *  
## pp.vac.1003          -0.111118   0.068877  -1.613 0.106682    
## pp.vac.1004          -0.242739   0.095917  -2.531 0.011383 *  
## school.close.red2     0.001774   0.043782   0.041 0.967688    
## school.close.red3     0.128579   0.036852   3.489 0.000485 ***
## pop_denlog           -0.084909   0.035125  -2.417 0.015635 *  
## restric.gather.red1  -0.223756   0.125291  -1.786 0.074117 .  
## restric.gather.red3  -0.189786   0.098053  -1.936 0.052924 .  
## restric.gather.red4  -0.235337   0.092128  -2.554 0.010636 *  
## time                  0.006673   0.005550   1.202 0.229177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef_model <- summary(model.glm)$coefficients$cond[,1][-1]

coef <- as.data.frame(coef_model)

# AME time
AME.time <- coef_model[12] * mean(dnorm(predict(model.glm, type="link")))
z12 <- amex(AME.time)


# pop.density 
AME.popden<- coef_model[8]* mean(dnorm(predict(model.glm, type="link")))
z8 <- amex(AME.popden)


# vac
vac0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,  
"pp.vac.100"="0",
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac1 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,  
"pp.vac.100"="1",
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac2 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,  
"pp.vac.100"="2",
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac3 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,  
"pp.vac.100"="3",
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac4 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,  
"pp.vac.100"="4",
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")

vac10<- vac1 - vac0
vac20<-vac2-vac0
vac30<-vac3-vac0
vac40<-vac4-vac0

AME.vac10 <- mean(vac10)
AME.vac20<-mean(vac20)
AME.vac30<-mean(vac30)
AME.vac40<-mean(vac40)

z2<-amex(AME.vac10)
z3<-amex(AME.vac20)
z4<-amex(AME.vac30)
z5<-amex(AME.vac40)
# 
# AME.vac <- coef_model[2]* mean(dnorm(predict(model.glm, type="link")))
# z2 <- amex(AME.vac)

# transport

trs0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"="0",  
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")

trs1 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"="1",  
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
trs10 <- trs1 - trs0
AME.trs10 <- mean(trs10)
z1 <- amex(AME.trs10)

#school close
sch0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,  
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"="0", 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
sch2 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"="2", 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
sch3 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"="3", 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= data0.anal1$restric.gather.red,  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")

sch20 <- sch2 - sch0
sch30 <- sch3 - sch0 
AME.sch20 <- mean(sch20)
AME.sch30 <- mean(sch30)
z6 <- amex(AME.sch20)
z7 <- amex(AME.sch30)

# restrict gather
gather0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= "0",  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")

gather1 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= "1",  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")

gather3 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= "3",  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")

gather4 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red, 
"pop_denlog"=data0.anal1$pop_denlog,  
"restric.gather.red"= "4",  
"test.policy"=data0.anal1$test.policy,  
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")

gather10 <- gather1 - gather0 
gather30 <- gather3 - gather0
gather40 <- gather4 - gather0

AME.gather10 <- mean(gather10)
AME.gather30 <- mean(gather30)
AME.gather40 <- mean(gather40)

z9  <- amex(AME.gather10)
z10 <- amex(AME.gather30)
z11 <- amex(AME.gather40)

# Test policy 

# test0 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red, 
# "pop_denlog"=data0.anal1$pop_denlog,  
# "restric.gather.red"= data0.anal1$restric.gather.red,  
# "test.policy"="0",  
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
# test1 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red, 
# "pop_denlog"=data0.anal1$pop_denlog,  
# "restric.gather.red"= data0.anal1$restric.gather.red,  
# "test.policy"="1",  
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
# test2 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red, 
# "pop_denlog"=data0.anal1$pop_denlog,  
# "restric.gather.red"= data0.anal1$restric.gather.red,  
# "test.policy"="2",  
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
# test3 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red, 
# "pop_denlog"=data0.anal1$pop_denlog,  
# "restric.gather.red"= data0.anal1$restric.gather.red,  
# "test.policy"="3",  
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
# 
# test10 <- test1- test0
# test20 <- test2-test0
# test30 <- test3-test0
# 
# AME.test10 <- mean(test10)
# AME.test20 <-mean(test20)
# AME.test30 <-mean(test30)
# 
# z12 <- amex(AME.test10)
# z13 <- amex(AME.test20)
# z14 <- amex(AME.test30)

z_merge <- rbind(z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11,z12)

term3 <- c(
"Closing public transport", 
"1-<5%",
"5-<10%",
"10-<30%",
"Vaccine coverage\n>=30%",
"School close at some levels",
"School close at all levels",
"Log of population density",
"Restrictions on gathering\n100-<1000",
"Restrictions on gathering\n10-<100",
"Restrictions on gathering\n<10",
"Time\n"
  )


p_load(broom.mixed)
sum_table <- tidy(model.mul, conf.int = T)

# create new table summary of model
new_table <- sum_table[2:13,]
 
new_table1 <- cbind(new_table, term3, z_merge)

new_table2 <- new_table1 %>% select(term, term3, AME, estimate, conf.low, conf.high)

new_table2$term3 <- factor(new_table2$term3, levels = new_table2$term3)


p.vac.period<- ggplot(new_table2, aes(x = term3, y = estimate,
                     ymin = conf.low, ymax = conf.high)) +
    geom_hline( yintercept = 0, color = 'Black' ) + ylim(-0.6,0.6) + 
    geom_linerange(size =1, color = 'steelblue') + geom_point(size=1.5, color ='steelblue') + coord_flip() + 
  ylab("Probit Transformation of Average Weekly Growth Rate ") + 
  xlab("Predictors") +
  theme_minimal() + 
  theme(strip.background = element_rect(NA), 
  axis.text = element_text(size = 20),
  axis.title = element_text(size = 20),
  legend.text = element_text(size = 10),
  legend.title = element_text(size = 10),
  strip.text = element_text(size = 8),
  legend.position = "top")
p.vac.period

ame.plot <- ggplot(data = new_table2, aes(x= term3, y = AME)) + geom_bar(stat="identity",fill="steelblue") + coord_flip() +
  geom_hline( yintercept = 0, color = 'Black' ) +
  geom_text(aes(label=round(AME,2)), hjust= 1.05, color="black", size=5)+
  ylab("Average marginal effect (%)") +
  ylim(-1,1) +
  xlab("") +
  theme_minimal() + 
  theme(strip.background = element_rect(NA), 
        axis.text.y=element_blank(),
  axis.text = element_text(size = 20),
  axis.title = element_text(size = 20),
  legend.text = element_text(size = 10),
  legend.title = element_text(size = 10),
  strip.text = element_text(size = 8),
  legend.position = "top")
ame.plot

cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig3_2021_vac_cate_Sep29.pdf",
          width = 20, height = 10, #inch
          onefile = TRUE, family = "Segoe UI")
merge_bar <- ggarrange(p.vac.period, ame.plot, ncol = 2, nrow = 1, widths = c(1.5,1))
merge_bar

# cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig3_2021_vaccine_Sep29.pdf",
#           width = 20, height = 10, #inch
#           onefile = TRUE, family = "Segoe UI")
# merge_bar <- ggarrange(p.vac.period, ame.plot, ncol = 2, nrow = 1, widths = c(1.5,1))
# merge_bar
dev.off()
## png 
##   2